obj=InputData();
obj.gMat=[1000 0.25];
obj.gNdt=[0,0;0,1;8,1;8,0];
obj.gElt=[4,4,1,4,3;2,4,1,2,3];
obj.gBctDisp=[1,1,1;2,1,1];
obj.gBctLoad=[2,0,6;3,0,2];

E=obj.gMat(1,1); u=obj.gMat(1,2);%弹性模量和泊松比
D0=E/(1-u.^2);
De =[D0,D0*u,0;D0*u,D0,0;0,0,D0*(1-u)/2];

Node=obj.getline(obj.gNdt);
gKA=zeros(2*Node,2*Node);
gFA=zeros(2*Node,1);
%计算并组集整体刚度矩阵
iLenElt=obj.getline(obj.gElt);
for i=1:iLenElt
    iEN = obj.gElt(i,2);
    viNodes=obj.gElt(i,3:iEN);% 获取单元的结点编号
    Ke=zeros(2*iEN,2*iEN);
    %计算单元弹性矩阵
    iMatID=obj.gElt(i,1);
%     %处理荷载边界条件
%     iBct=obj.getline(obj.gBctLoad);
%     for ib=1:iBct
%         iNode = obj.gBctLoad(ib,1);%编号
%         for i=1:2
%             irow=(iNode-1)*2+i;
%             gFA(irow) = gFA(irow)+gBctLoad(ib,i+1);
%         end
%     end
    
    iBct=obj.getline(obj.gBctLoad);
    for i=1:iBct
        %获取节点坐标
        xNode=gBcload1(i,1:2);
        yNode=obj.gNdt(xNode(:),:);
        %计算长度
        length=sqrt((yNode(1,1)-yNode(2,1))^2+(yNode(1,2)-yNode(2,2))^2);
        %计算等效节点荷载
        if(gBcload1(i,3)==0)
            for j=1:2
                gFA(xNode(j)*2-1)=gFA(xNode(j)*2-1)+0.5*length*gBcload1(i,4);
            end
            sub=gBcload1(i,4)-gBcload1(i,5);
            gFA(xNode(1)*2-1)=gFA(xNode(1)*2-1)+1/3*length*(-sub);
            gFA(xNode(2)*2-1)=gFA(xNode(2)*2-1)+2/3*length*(-sub);
        else
            for j=1:2
                gFA(xNode(j)*2)=gFA(xNode(j)*2)+0.5*length*gBcload1(i,4);
            end
            sub=gBcload1(i,4)-gBcload1(i,5);
            gFA(xNode(1)*2)=gFA(xNode(1)*2)+1/3*length*(-sub);
            gFA(xNode(2)*2)=gFA(xNode(2)*2)+2/3*length*(-sub);
        end
    end
    
    %处理位移边界条件
    iBct=obj.getline(obj.gBctDisp);
    for i=1:iBct
        iNode = obj.gBctDisp(i,1);
        for j=1:2
            irow=(iNode-1)*2+j;
            if (obj.gBctDisp(i,j+1)==1)
                gKA(irow,:)=0.0;
                gKA(:,irow)=0.0;
                gKA(irow,irow)=1.0; %划零置一法
                gFA(irow)=0.0;
            end
        end
    end
    %求解方程
    vdDisp=gKA\gFA;
    %记录结点位移
    gNTu=zeros(iNdt,2);
    for id=1:iNdt
        gNTu(id,1:2)=vdDisp(2*id-1:2*id);
    end
    %显示结点的计算位移值
    disp('nodal displacement:x, y');
    disp(gNTu);
    %**************************************
    %在屏幕上绘制变形图
    hold on %绘图时不清屏
    dFau=0.2; %变形图放大系数（调整此系数显示合理变形图）
    %绘制变形前网格
    trisurf(objgElt(:,3:5),obj.gNdt(:,1),obj.gNdt(:,2),zeros(size(gNdt,1),1))
    view(2); axis equal; axis off; axis tight;
    pause(2.0); %暂停
    %获取变形后数据
    DDisp=obj.gNdt+gNTu*dFau;
    trisurf(obj.gElt(:,3:5),DDisp(:,1),DDisp(:,2),zeros(size(DDisp,1),1));
    view(2); axis equal; axis off; axis tight;
    clear;
    %计算单元应变矩阵 对三角形单元
    %获取单元的结点坐标3*2矩阵
    xe=obj.gNdt(viNodes(:),:);
    %计算形函数系数（ ai 、 bi 、 ci
    for i=1:3
        j = i+1;
        if(j>3)
            j=j-3;
        end
        k = i+2;
        if(k>3)
            k=k-3;
        end
        ai(i) = xe(j,1)*xe(k,2)-xe(k,1)*xe(j,2); %ai/aj/am
        bi(i) = xe(j,2)-xe(k,2); %bi bj bm
        ci(i) = xe(k,1)-xe(j,1); %ci cj cm
    end
    delta=sum(ai); % 单元面积 *2
    dA=delta/2;
    if(delta<=1e-10)
        error(' 单元面积小于 0!');
    end
    % 计算应变矩阵 Bm
    Bm=[bi(1) 0 bi(2) 0 bi(3) 0;
            0 ci(1) 0 ci(2) 0 ci(3);
            ci(1) bi(1) ci(2) bi(2) ci(3) bi(3)];
    Bm=Bm/delta;
    % 计算单元刚度矩阵
    Ke=Bm'*De*Bm*dA;
    % 组集单元刚度矩阵 Ke 到整体刚度矩阵 gKA
    iLR = 2*iEN; % 单元的总自由度数
    viDof = zeros(2*iEN,1); %各结点的总自由度编号
    for il=1:iEN
        inod = viNodes(il);
        viDof(2*il-1:2*il) = 2*inod-1:2*inod;
    end %for il
    % 组集
    for ir=1:iLR
        irg=viDof(ir);
        for ic=1:iLR
            icg=viDof(ic);
            %ir/ic 为单元刚度矩阵的局部编号
            %irg/icg 为其对应的总体刚度矩阵编号
            gKA(irg,icg)= gKA(irg,icg)+Ke(ir,ic);
        end
    end
end